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Abstract The status of the hydrodynamical modelling of nonlinear multi-mode 
stellar pulsations is discussed. The hydrodynamical modelling of steady 
double-mode (DM) pulsations has been a long-standing quest that is 
finally being concluded. Recent progress has been made thanks to the 
introduction of turbulent convection in the numerical hydrodynamical 
codes which provide detailed results for individual models. An overview 
of the modal selection problem in the HR diagram can be obtained in 
the form of bifurcation diagrams with the help of simple nonresonant 
amplitude equations that capture the DM phenomenon. 
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1. INTRODUCTION 

Self-excited multi-mode pulsations are quite common among the 
less evolved luminous stars, such as the delta Scuti stars, but they 
mostly involve nonradial modes of oscillation. In contrast, the clas- 
sical variable stars, viz. the Cepheids and RR Lyrae stars, are be- 
lieved to be radially pulsating, and two independent frequencies are 
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typically identified in the Fourier spectrum. These pulsating stars 
are referred to as double-mode even though, in principle, there could 
be more than two modes involved with locked frequencies. (We 
should add that there is evidence that suggests the presence for 
nonradial modes as well (Kovacs et al. , Moskalik, Olech et al. in 
Szabados & Kurtz 2000) . 

In the Galaxy, the Beat Cepheids, also called double-mode (DM) 
Cepheids are relatively rare. Only about a good dozen are known. 
On the other hand, the recent observations of the Magellanic Clouds 
(Beaulieu et al. 1995, 1997, Welch et al. 1995, Udalski et al. 1999) 
have shown that in these galaxies Beat Cepheids are quite com- 
mon. Beat Cepheid pulsations occur either in the fundamental/first 
overtone (F/Oi) modes, or in the first/second overtone (O1/O2) 
modes. In the Galaxy there is only one O1/O2 known DM pulsator 
(CO Aur), but in the SMC Udalski et al. 1999 claim 70 candidates 
out of 93 DM Cepheids. 

The occurrence of DM RR Lyrae stars (RRd stars), depends on 
the type of cluster. For example, not a single RRd star is observed 
in 00 Centauri (cf. Kovacs in this Volume) while in M68, 9 RRd 
stars are known out of 37 RR Lyrae stars. 

The DM pulsators are important for pulsation theory in that 
they impose very stringent requirements on the numerical mod- 
elling. Double periodicity and indeed the overall view of the modal 
selection picture provide many more observational constraints than 
their single-mode siblings. 

Despite the frequent occurrence of DM pulsations in nature, the 
numerical modelling of this type of pulsation had remained a seri- 
ous challenge until recently. In fact, it had become abundantly clear 
that purely radiative models, i.e. models that disregarded connec- 
tive transport, were not capable of yielding DM pulsations, except 
on a purely transient basis, i.e. they were switching from one mode 
to another, but too fast to account for the observed fairly steady 
nature of DM pulsators. Actually, some purely radiative models of 
RRd stars had been found already by Kovacs and Buchler (1993), 
but those results were not satisfactory. However, recently, DM be- 
havior has been found simultaneously, and fully independently, in 
RR Lyrae models by Feuchtinger (1998) and in Cepheid models by 
Kollath et al. (1998), this with different numerical methods, viz. the 
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Vienna and the Florida codes. The breakthrough came as a result of 
the inclusion of time-dependent turbulent convection in the models. 

2. AMPLITUDE EQUATIONS - SIMPLE 
MODELS FOR DM PULSATIONS 

Numerical hydrodynamics may provide us with an accurate de- 
scription of the pulsations of an individual stellar model, but it does 
not give us an overview of the behavior when model parameters are 
changed. On the sole basis of numerical hydrodynamical modelling 
one often has no clue why neighboring models may exhibit very 
different pulsational behavior. Examples of such sensitive behavior 
will be discussed below. 

Therefore a global way had to be devised for gaining a view of 
the global modal selection picture (or the bifurcation diagram as it 
would be called more generally). The simplest global way to de- 
scribe the interaction of pulsation modes, including DM behavior, 
is through the use of amplitude equations, which are also known 
as normal forms. A derivation of the amplitude equation formal- 
ism applied to stellar pulsation can be found elsewhere (Buchler 
& Goupil 1984, Buchler & Kovacs 1986, cf. also Dziembowski & 
Kovacs (1984), Takeuti (1985), Goupil & Buchler 1994 and for a 
review Buchler 1993). This formalism is based on dynamical sys- 
tems theory, and it is very general and fundamental. The ampli- 
tude equations give the temporal behavior of the amplitudes and 
the phases of the modes that are involved in the nonlinear pulsa- 
tion. They are of fundamental importance for understanding the 
changes in the morphology of the Fourier decomposition coefficients 
of the light and radial velocity curves as a function of pulsation 
period. For example, they explain the nature of the well-known 
Hertzsprung progression of the classical Cepheids. The amplitude 
equations show unambiguously how these phenomena are related 
to and caused by internal resonances of the pulsating mode with an 
overtone (for a review cf. Buchler 1993). 

We consider here the situation in which only two modes are in- 
volved in the pulsation and there exists no resonance condition be- 
tween them nor with the other modes of the star. This is in contrast 
to earlier purely radial modelling which had required the presence 
of resonances (Kovacs & Buchler 1993), but which turned out to be 
unsatisfactory. 
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We denote the linear eigenvalues by <jj, for an assumed exp(cxt) 
dependence, 
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The absence of resonance means that there is no relation of the 
form riiUo + n^x ~ or niu + n 2 Ui + n 2 cJ3 0, where rii, 112, 
are small positive or negative integers. 

The generic amplitude equations, appropriate to the nonresonant 
situation, when truncated at 5th order read (e.g. Buchler 1993): 

d = a (a + Qoo\ao\ 2 + <5oi|ai| 2 + So|a | 2 |ai| 2 

+ J R ok| 4 + #oih| 4 ) (1.4) 

di = ai(<7i + Qio|a | 2 + <5n|ai| 2 + S , i|a | 2 |ai| 2 

+ J R 10 |a | 4 + J Rn|a 1 | 4 ) (1.5) 

The a/s are the complex amplitudes of the two excited modes, and 
Qjk, Sj and Rjk are the complex nonlinear coupling constants. 

Note that only selected powers of the amplitudes appear in Eqs. 
(1.4,1.5); for example there are no quadratic terms. Normal form 
theory (e.g. Coullet & Spiegel 1983) shows (a) that these equa- 
tions are generic, i.e. they apply to any system that is character- 
ized by having two excited nonresonant modes, be it a pulsating 
star, a biological system, or any other system, in which the relative 
growth-rates of the modes are small (rjj <^ 1); and (b) that the 
omitted terms are not essential for the modal selection (bifurcation 
diagram), i.e. that the nature of the dynamical behavior does not 
depend on them. 

Usually it is more convenient to use real amplitudes Aj(t) and 
phases (pj(t), defined by a,j(t) = Aj(t) expiipj(t). Inserting this 
relation into the above equations one gets: 

A i = A 3 + qMl + qjl A\ + SjAqA\ + r j0 Al + rjl A\) (1.6) 
and 

<Pj = + QjoAl + m A l + h A l A l + r j0 A% + f^Aj, (1.7) 
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where we have introduced the real constants by the relations: 

Qjk = Qjk + iqjk (1-8) 
Sj = Sj + isj (1.9) 
Rjk = r jk + if jk (1.10) 

One notes that the introduction of the real amplitudes Aj in Eqs. (1.4, 
1.5) produces a complete decoupling of the amplitudes from the 
phases. 

In principle, the constants Q, S and R can be derived from the 
stellar structure equations and the linear eigenvectors as was shown 
in Buchler & Goupil (1984). In practice, however this is a daunt- 
ing task that has only been attempted in some very specific cases 
(Takeuti & Aikawa 1981), Klapp et al. 1985, Dziembowski & Kro- 
likowska 1985). An alternative approach that works quite well is to 
derive them from numerical hydrodynamical studies (e.g. Buchler 
& Kovacs 1987, and below). 

In the earliest use of amplitude equations for describing DM 
pulsations (Buchler & Kovacs 1986) the amplitude equations were 
truncated at the lowest nontrivial, i.e. cubic, terms. It was shown 
that in this order the single-mode fixed points and steady DM pul- 
sations cannot simultaneously be stable for the same stellar model. 
Subsequently, the behavior of the hydrodynamical models and the 
observational constraints imposed by the Beat Cepheids and RR 
Lyrae stars have forced us to consider also the next order, quintic 
terms (Buchler, Yecko, Kollath & Goupil 1999, hereafter BYKG). 
The properties of these equations and their applicability to the DM 
stellar pulsators were discussed in that paper. Here we present a 
slightly different discussion and approximation. While, for simplic- 
ity, BYKG retained only the r 00 and r n cubic terms, instead, we 
keep the Sq and s± cubic terms because they are sufficient to give a 
good, and in fact a better fit to the results of the hydrodynamical 
calculations. 

Time-Independent Amplitude Equations 

If our primary interest is the study of steady pulsations, i.e. oscil- 
lations with constant amplitudes, then we need to consider only the 
fixed points (FP's) of the amplitude equations Eq. (1.6), defined by 
A = Ai = 0. When only one FP amplitude is nonzero (a single- 
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mode FP) the corresponding full-amplitude (nonlinear) pulsations 
of the star have constant amplitudes and they are mono-periodic. 
These pulsations are then called limit-cycles, when in addition they 
are stable. 

Eqs. (1.7) show that for FP's the phases become linearly increas- 
ing functions of time (because the A^s are constant) and they thus 
yield the nonlinear modal frequencies cuj = <pj. 

The possible FP's are therefore obtained by looking for the solu- 
tions of the 4 combinations of the equation pairs: 

K + q 00 A 2 + q Q1 Al + s A 2 Al + r m Al + r Q1 A\ = (1.11) 

A 1 = (1.12) 

Ki + <?io^o + 9n^i + SiA 2 Al + r 10 A* + r lx A\ = (1.13) 

A = (1.14) 

Without loss of generality, we label the two modes and 1, because 
most of the time we will be talking about the fundamental and first 
overtone modes. 

Eqs. (1.11 - 1.14) are a function of the squared amplitudes. It is 
therefore sometimes advantageous to work in (Aq, A\) space. Thus 
in Figure 1.1 we have plotted the loci defined by the first of each of 
Eqs. (1.11, 1.13) as solid and dotted lines, respectively, for an RR 
Lyrae model that will be discussed below. The remaining two loci 
are the positive x and y axes. The solid and open circles represent 
stable and unstable FPs which are the intersections of these loci. 
Note that the situation is much simpler when all the quintic terms 
are disregarded as in Buchler & Kovacs (1986) because then the 
loci degenerate into a set of straight lines, and clearly there can be 
at most one DM solution in that case. Here we consider the more 
general case. 

The trivial FP of the amplitude equations is the static stellar 
model, i.e. with both A = and A 1 = 0. If both growth-rates are 
negative (/c , «i < 0) the FP is stable, and so is the star. The next 
simplest case is where only one k is positive in which case the only 
possible pulsation state of the star is a limit-cycle in that mode. 
For the existence of a DM solution both growth-rates need to be 
positive, but as we will see it is not a sufficient criterion for steady 



Double-Mode Stellar Pulsations 7 



Table 1.1 Possible fixed point structures of the nonresonant amplitude equations (in 
addition to the unstable trivial FP, Ao = A\ = 0), assuming ko, «i > 0. S: stable, U: 
unstable. 



Number of FP's 


F 


Oi 


DM1 


DM2 


2 


S 


u 






2 


u 


s 






3 


s 


s 


u 




3 


u 


u 


s 




4 


s 


u 


s 


U 


4 


u 


s 


s 


u 



DM oscillations. With the nonresonant amplitude equations only 
a limited number of possible limit-cycle combinations are possible, 
and they are listed in Table 1.1. 




Figure 1.1 DM solutions of the amplitude equations in the (Aq, A\ plane for an RR 
Lyrae model. Solid circles are stable and open circles are unstable fixed-points 
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Figure 1.2 DM solutions of the amplitude equations, see text Eq. (1.19); Solid circles 
denote stable and open circles denote unstable DM solutions. 

In order to make the problem tractable we now keep the s terms, 
but omit the r quintic terms in Eq. (1.11,1.13). It is simple to find 
the fundamental and first overtone limit-cycle amplitudes 



A 2 



-Kq 



,lc 



<?00 

9n 1 



A 1M = 
A i c = 0. 



(1.15) 
(1.16) 



Inserting those values into the amplitude equations corresponding 
to the other mode, we can easily obtain the conditions for the linear 
stability of the fundamental and overtone limit-cycles, respectively: 



«! = «! + q w A 2 1c = ki- Kq— < 

<?oo 

- — i A 2 901 ^ n 

K = K + goiA ic = «o - «i — < 

<?11 



(1.17) 
(1.18) 



We illustrate the modal behavior with an RR Lyrae model for 
which we have computed the nonlinear coupling coefficients (see 
Table 1.2). The numbers have been rounded and the g's and s's 
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have been normalized to set the limit-cycle amplitudes to the unity 
(e.g. g o = -«o and q u = -Ki). 

In order to see the behavior of the solutions along a sequence of 
RR Lyrae models, we fix all these parameters but let k$ vary. This 
sequence of solutions, obtained from Eqs. (1.11 - 1.14), has already 
been depicted in Fig. 1.1. 

From the stability of the limit-cycles one further finds that the 
range of stable fundamental mode pulsations is k,q > 0.007. The 
range for the stable overtone is k < 0.01. This indicates that 
both limit-cycles are possible (the 'either-or' regime in the pulsation 
jargon) and these multiple solutions exist in the 0.007 < k < 0.01 
interval, and necessarily there is a DM fixed point for that regime 
which is furthermore unstable. 

In Figure 1.1 one easily misses some important and interesting 
solutions that occur in an extremely narrow range of k = 0.0100 
to 0.0101, located between the last two subfigures. In fact for the 
parameters of Table 1.3 the upper, solid curve undergoes two inter- 
sections with the dotted curve which corresponds to two additional 
DMs. 

In order to exhibit this behavior in a more apparent way we 
consider the quantities Fj, derived from Eqs. (1.11, 1.13). 

Plots of 5F = Fi — Fq as a function of A$ give the amplitudes of 
the DMs solution. It can be seen than the fixed point is stable if the 
slope of F is steeper than that of i.e. ddF/dA^ > 0. Otherwise 
no stationary DM pulsation is possible. 

A stable DM solution exists and both limit-cycles lose their sta- 
bility if and only if 5F is an increasing function of Aq in the neigh- 
borhood of the fixed point. With the opposite slope both the fun- 
damental and the first overtone limit-cycles are stable, giving the 
"either-or" region of pulsation (Buchler & Kovacs 1986). 

Figure 1.2 shows how the SF = F\ — F Q function varies as k 
is changed. Only a very narrow range exists where the nearly 
parabolic function has two intersections with the zero line. 

Both Figs. 1.1 and 1.2 show that the DM amplitudes are al- 
ways smaller than the corresponding SM limit-cycle amplitudes. 
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Table 1.2 Growth- rates and nonlinear coupling coefficients from an RR Lyrae model. 
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Ki 




qio 
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0.042 


-0.01 
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-0.042 
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-0.001 



When the DM occurs near the A\ axis, then A\^ m < Aij c , but 
Ao y dm *C Aqj c . With the given parameters the maximum funda- 
mental amplitude is only 30% of the limit-cycle amplitude. For 
most of the RRd stars the fundamental mode amplitude is indeed 
smaller than that of the first overtone (see Kovacs in this Volume) , 
indicating that the above RR Lyrae model is on the right track. 

Time-Dependent Amplitude Equations 

The time-dependent solutions of the amplitude equations become 
essential when we interpret the numerical hydrodynamical integra- 
tions or explain the changes due to e.g. evolutionary changes of the 
star. As we see above, very small changes in /c's can induce signifi- 
cant changes in the possible pulsation state of the star. In our test 
problem if the fundamental mode growth-rate is slowly increasing, 
at k,q = 0.01 the overtone mode loses its stability. Then the model 
evolves smoothly to the DM state, if the system has been pulsat- 
ing in the overtone mode before. With a further increase of k,q, at 
k 0.01007, the fundamental mode pulsation remains the only 
stable state. Then the amplitude equations can be used to esti- 
mate how the stars switches from the DM state to the fundamental 
limit-cycle. 

In Fig. 1.3 the temporal variation of the amplitudes are shown 
during the simulated mode switching. The integration of the ampli- 
tude equations was initiated (at t=0) at the stable DM fixed point, 
and Kq was continuously increased at a 10~ 7 years -1 rate. The DM 
oscillations became unstable during the first 10 years, but it took 
more than 150 years to evolve the star to the fundamental mode 
pulsations. 

Right at the bifurcation point where the stable DM pulsation 
trades stability with the unstable fundamental SM limit-cycle the 
growth-rate of the DM pulsation is infinitely small and it takes 
evolution for the star to move into the regime where this growth-rate 
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achieves sufficiently large values for the switch to the fundamental 
limit-cycle to occur. It is therefore clear that the mode switching 
time-scale depends both on the evolution time-scale and on the 
sensitivity of the modal growth-rate to the evolution time-scale. 
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Figure 1.3 Simulated mode switching as predicted by the amplitude equations, see 
text; solid: A , dotted: Ai. 



In section 4. we discuss how the hydrodynamical calculations can 
be used to reconstruct the amplitude equations. 

3. HYDRODYNAMICAL MODELLING OF 
DOUBLE-MODE PULSATION 

The Hydrodynamical Code 

The classical variable stars that we consider in this review are 
radial pulsators. The global motions thus can be modelled by 
ID hydrodynamics (see however Kovacs in this Volume for indi- 
cation of nonradial modes in RR Lyrae stars). Since real 3D mod- 
elling remains a dream in nonlinear pulsation calculations, the local 
nonradial flows i.e. turbulence have to be treated in some one di- 
mensional approximation. Recipes and pulsation codes for ID tur- 
bulent convection (TC) have been developed by different authors 
(e.g. Stellingwerf 1982, Kuhfufi 1986, Gehmeyr & Winkler 1992, 
Feuchtinger, M. U. 1998, Yecko, Kollath, & Buchler 1998). 

The fluid dynamics part of the model calculations are given by 
the following equations: 

du Id, s GM r , 

~m = —jr(p+ Pt+Pv) — 1-20 

at p or r 2 
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de dv Id 



dt dt pr- 2 Or 



,2 



(F r + F c ) + C (1.21) 



The turbulent motion of the gas and the convection interacts 
with the hydrodynamics of the radial motion through the convec- 
tive flux F c , the viscous eddy pressure p v the turbulent pressure p t , 
and finally, through an energy coupling term C. Those terms can be 
derived by averaging the 3D hydrodynamics equations. This proce- 
dure introduces a sequence of moment equations (see e.g. Canuto 
1998; Canuto & Dubikov 1998). In the simplest recipe we truncate 
all this equations, but leave the time-dependent diffusion equation 
for the turbulent energy e t : 

de t , s d v 1 d ( , „ \ 

- + (pt+Pv) - = ---(SF,)-C (1.22) 

The coupling term that connects the gas and the turbulent energy 
equations is given by: 

1/2 

C = -^-« d (e t - S t ) , (1.23) 

where A = a A H p) H p = pr 2 / (pGM) is the pressure scale height «a 
is the mixing length parameter and ad is a dimensionless parameter. 

Both the convective flux and the source term of the turbulent 
energy (S t ) depend on the entropy gradient: 

H v ds , 

Y = v - — (1.24) 

c p or 

S t = (a s a A ) 2 ^[3TYf pec , (1.25) 



F c = a c a A pe] /2 c p T Y f pec , (1.26) 

where (3 is the thermal expansion coefficient a s , a c are parameters 
and other symbols have their usual meanings. The Peclet correction 
f pec accounts for the decrease of convective efficiency when radiative 
losses are important. We approximate this correction factor by 
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(1.27) 



pec 



1 + a r Pe~ l 



D r 



Pe 




(1.29) 



(1.28) 



(1.30) 



The remaining quantities 



are defined as 



Pt = 



a p pe t 



(1.31) 



Vv = 



Ft = 




(1.32) 



(1.33) 



The TC model equations are based on a very simple physical 
picture of turbulent energy generation and diffusion, and convective 
energy transport. The various expressions can essentially be derived 
from dimensional arguments (even though numerous attempts have 
been made to derive them from more detailed analyses - for a dated, 
but still excellent review, see Baker 1987). The above definitions 
of the turbulent energy source and of the convective flux are not 
unique, and slightly different expressions have been used by different 
authors. The formulation adopted is of that of Gehmeyr & Winkler 
(1992). For a comparison of the different recipes used for stellar 
pulsation calculations see Buchler et al. (1999). 

The dimensionless a parameters of our turbulent convection model 
equations are of order unity, but theory provides no guidance as to 
their numerical values. We will resort to a comparison of numerical 
results with empirical data, such as the position and the width of 
the instability strip, to calibrate those parameters. 

We have not yet performed a full calibration, which is a giant task 
because of the number of parameters, and in fact it is not even clear 
whether all observational constraints can be fulfilled with the use of 
such a simple model for turbulent convection. In the meantime we 
have chosen several sets of a parameters that satisfy some of the 
constraints, and we have explored the consequences of turbulent 
convection for the linear and nonlinear properties of pulsation. 
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It has been comforting that, for the first time, we have been able 
to calculate DM pulsation with these parameters for a wide set of 
(physical) stellar parameters. However the results still miss some of 
the observational constraints, e.g. the model amplitudes are small 
compared to the observed ones, the modelling of 1st overtone/2nd 
overtone DM pulsations is still not quite satisfactory. At this time 
no detailed or comprehensive comparison has been be made between 
the observations and the DM modelling. Thus this paper is a status 
report on modelling rather than a closed and self contained theory 
of DM pulsations. However, the tendencies, e.g. along a sequence 
with different temperatures, can be used as primary tests for our 
models. 

The first step in the numerical modelling is the construction of the 
static envelope model, and its linear stability analysis which yields 
the possible pulsation periods and the corresponding excitation- 
rates. These modal growth-rates tell us only the rate of increase 
or decrease of the oscillation for tiny amplitudes. Furthermore, if 
more than one mode is excited we cannot infer from this information 
which of the modes will grow into full-amplitude nonlinear oscilla- 
tions. This 'modal selection' problem has already been discussed 
in section 2. where it was shown how it can be attacked with the 
amplitude equation formalism. 

There exists no really simple way for computing the limit-cycles 
or the final DM or multi-mode pulsation state of the stellar mod- 
els. The hydro dynamical modelling is an initial value problem, 
and we need to evolve the model directly with the hydrodynamics 
code, or through a relaxation method to obtain the stable oscil- 
lations. Furthermore, the model can evolve to different final sta- 
tionary states depending on the initial kick of the static envelope. 
Fortunately, the investigated models have only one or two possible 
stable pulsation states, since for radially pulsating envelopes only a 
few modes are generally linearly unstable (at least for the Cepheids 
and RR Lyrae). In the case of nonradial pulsations, where dozens 
of modes are linearly unstable, the number of stable NL oscillations 
can be higher. 

A great advantage of the Florida code is that it can relax to the 
periodic solutions by iteration and it provides the stability (Flo- 
quet) analysis of the limit-cycles. When both the fundamental 
mode and the first overtone are linearly unstable but the other 
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modes are stable, then the Floquet exponents of their limit-cycles 
provide the primary information on modal selection. 

DM Solutions in the Numerical Calculations 

The computation of the possible full-amplitude (or nonlinear) pe- 
riodic pulsation states of a stellar model is essentially routine work 
with the relaxation code. However there exists no similar method 
for finding DM or other multi-mode full amplitude pulsations. Even 
when the fundamental and the first overtones are both linearly un- 
stable, nothing guarantees the existence of steady DM pulsations. If 
both periodic limit-cycles are unstable against perturbations, then 
the model should oscillate in a DM state only. However, as the 
solutions of the amplitude equations indicate (see Section 2.) stable 
DM solutions can exist even when one of the single-mode limit- 
cycles is stable. In such a case the Floquet analysis is insufficient 
for identifying DM pulsations. 

The most efficient way of kicking the initial static model with the 
velocity eigenvectors from the linear stability analysis. In order to 
obtain different transient paths we thus use linear combinations of 
the fundamental and the overtone mode eigenvectors with a surface 
amplitude of rs 10 km/s. 
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Figure 1.4 The initial transient in RR Lyrae model pulsations 

In Figure 1.4 the radius variation of an RR Lyrae model is shown 
for the initial 100 days of model integration. The multi-component 
nature of the signal can be seen from an inspection of the curve 
- the beats are clearly visible. The amplitudes of the modes are 
seemingly saturated, but further integration of the model shows 
that this is not the case. At that stage of the hydrodynamical 
evolution one cannot predict the final pulsation state: it can be a 
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limit-cycle with any of the unstable modes, or a DM pulsation. The 
system generally has to be integrated for a very long time (usually 
for thousands of cycles) to determine whether the amplitude of the 
modal components saturates at some finite value. 

To emphasize the importance of long term integration we have 
calculated the amplitudes from an extended continuation of the 
transient shown in Fig. 1.4. The variation of the amplitudes is dis- 
played in Fig. 1.5. Note the different time-scales on the two figures: 
the curve shown in Fig. 1.4 corresponds to the very beginning of 
the amplitude evolution on Fig. 1.5. 




40 60 
time [years] 



Figure 1.5 The evolution of the amplitudes after the transient shown in Fig. 1.4. 
Solid line: fundamental mode, dotted line: first overtone. 



Claiming DM behavior is therefore not an easy task. Very long 
transients with mixed-mode behavior exist in the hydrodynamical 
calculations, even when the model finally converges to a SM limit- 
cycle. This means that the knowledge of the Fourier spectra and 
of the velocity and light-curves from a single hydrodynamical cal- 
culation of a model cannot provide a sufficient proof of steady DM 
behavior. A more elaborate set of calculations is necessary. This is 
the topic of the next sections. 

Time-Dependent Amplitudes and Phases 

Frequently a visual inspection of the time-series that is obtained 
from the hydrodynamical calculations may indicate that the pulsa- 
tions have reached a steady state. However these indications can be 
quite treacherous. Even the Fourier analysis of the light, velocity or 
radius variation can be misleading because of the temporal changes 
of amplitudes and periods in the transitory stage of model calcu- 
lations. In order to obtain reliable information from the numerical 
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results it is necessary to apply a well suited time-frequency analysis 
and derive instantaneous amplitudes and periods. 

A variety of such methods are available. For example, Kovacs, 
Buchler & Davis (1987) used a time-dependent Fourier method. 
They performed linear least-squares fits with sine functions to very 
small successive portions of the data. The bases of the successive 
fits were shifted by small portions of the fit, and the results were 
averaged to get smooth results. The disadvantage of this method is 
that short bases provide large errors on the amplitudes, while longer 
bases introduce temporal averaging destroying the instantaneous 
nature of the amplitudes and phases. 

Instead of the time-dependent Fourier analysis one can follow 
Gabor (1946) to reconstruct the time- dependent amplitudes and 
frequencies of the pulsation modes. Let s(t) represent the real part 
of an assumed complex analytical function a(t). The imaginary 
part s(t) of a(t) can then be obtained via a Cauchy integral, which 
through contour deformation becomes a Hilbert transform. Physi- 
cists are generally familiar with this analytic signal concept through 
the Kramers-Kronig dispersion relations (Jackson 1975). 



It is easy to verify for example that with s(t) = Acosut one finds 
s(t) = A smut. The Gabor construction makes it possible to unam- 
biguously define the phase tp(t) and the amplitude A(t) of a signal. 
Consequently it allows one to extract the instantaneous frequency 
(Cohen 1994). 

It is possible to extend the method to multi-component signals. 
To obtain the instantaneous amplitudes and periods of the pulsation 
modes separately one first filters the signal to eliminate the power 
from the other frequency components. It is convenient to make the 
filtering in Fourier space, and combine it with the definition of the 
analytic signal 




(1.35) 



(1.34) 



A(t)e Mt) 



(1.36) 
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1 roo rco 

Z k (t) = a k {t)e 1 ^ = - / H(u - oo k ) / s(t')e M '^ Ut'du, 

IX JO J-oo 

(1.37) 

where H{uj — u k ) is the window of the filtering, centered on u k . A 
Gaussian window with a half width of 0.2 c/ d gives a satisfactory re- 
sult for RR Lyrae models. The resulting amplitude and phase gives 
the temporal evolution of the given mode. The least-squares fitting 
of a Fourier sum to the radius variation gives essentially the same 
results. However the analytic signal gives better resolution and it 
is not necessary to average the amplitudes in time to get smooth 
results. Thus it is more appropriate for further analysis. A great 
advantage of the analytic signal is that the resulting amplitudes can 
be directly used to fit the coefficients of the amplitude equations 
in a differential form - a simple linear least-squares method. The 
amplitudes calculated by the method of KBD87 are not suitable for 
a differential fit, but only for a complicated nonlinear least-squares 
fit to the integral curves. 

Modal Selection Problem 

Armed with this tool we can now analyze the modal selection 
problem in the hydrodynamical models. A time consuming, but 
very reliable method is to perform a number of hydrodynamical 
integrations with different initial perturbations of the same static 
model, and then to extract the slowly varying amplitudes with the 
help of a time-frequency analysis. The resulting phase-portrait 
(Ai(t) vs. Ao(t)) gives a clear global overview of the nonlinear 
behavior of the model, and thus permits one to see immediately if 
stable DMs exist. 

In Fig. 1.6 we show a representative set of such hydrodynamical 
tracks for a Cepheid model. The initial perturbations have been 
consisted of a mixture of the fundamental and the first overtone 
eigenvectors. Fig. 1.6 provides conclusive proof of the presence of 
stable beat pulsation with amplitudes A = 0.016 and A 1 = 0.028 
that coexists with the stable fundamental mode with amplitude 
Aq = 0.049. In this case there is hysteresis and the actual state 
of the star depends on its former evolutionary path. From Fig. 1.6 
one can also infer that two additional nontrivial fixed points exist 
for the given model in the phase-space of the amplitudes, namely a 



Double-Mode Stellar Pulsations 19 



0.03 



0.02 



<0 



0.01 





I 



if 



i 

f'» 

J J w 

J I -i J 

I ]•/ I 

'I I.i /' 

i /.I I 

i i i r 

i r i f 

lilt 
I J 1 V 
/./ /./ 

r i f i 
till 
•i n i 
•i 1'J /. 

1 J J A 
I.I /•/ 
£ / /. / 
.7 /.*/ / 

/V / ^ 



///////// 
J 1 1 I I ! I I I I 
////////// 
I 11 I I I 1 ! 1 
11111111 
! I I 1 I I 
kjtlll 

V / / / 
. . i IV// 
i:i i i iT^ ^ 
.,/./// ill \^ 
i i i'i i i ill i 
1 l f l l l i l l 
1 1 J 1 1 I j 1 1 I 
I I J I I 1'J J 1 1 
1 i'i 1 i ill 1 i I 
i C i i i + i / i I 
i 7 1 //.•// 1 1 1 
i.i i i /. / / i g/ 

/.///¥// /// / 
////'/// /• / / 
J // S ////// / 
///'/////// 
/////£/ /v/ 
/ X'/ / /,•/ / fy /, 
/.V / y.'y y ^y y 
'y y &'y y S / y 




0.01 



0.02 



0.03 



0.04 



0.05 



A [(5R/RJ 



Figure 1.6 The flow in the (Ao, Ai) phase-space. The large dots represent the hy- 
drodynamical results. The short lines represents the flow in the phase-space. The 
two separatrices are clearly visible: the first runs from the origin upward between the 
fourth and fifth evolutionary tracks, and the second one is the arc connecting the two 
single-mode fixed points. 



first overtone limit-cycle and a DM saddle-point, but these states 
must be unstable. The observed behavior of the model is in a good 
agreement with the results found with the amplitude equations. 
The amplitude equations can be used to find all the fixed points, 
even if the hydrodynamical trackss do not actually reach them, 
either because of insufficient integration time, or because they are 
unstable and thus simply unreachable this way. The flow field that 
is obtained from the amplitude equations is represented by the short 
thin lines in Fig. 1.6. 

The procedure that we have described, namely numerical hy- 
drodynamics integrations with judiciously chosen different initial 
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conditions, combined with a time-frequency analysis and with a 
Floquet analysis provide a complete picture of the possible pulsa- 
tion states of a given model and of their stability. By repeating the 
calculations with different effective temperatures and luminosities 
one can thus map the whole picture of modal selection on the HR 
diagram. 

4. RESULTS 

We have already mentioned that the calibration of the turbulent 
parameters is a big task that has not been finished yet. The results 
that we present here can show the general tendencies that exist in 
the models, but we cannot pretend that the models can be used to 
infer for example physical parameters from individual variable star 
observations. 

The Coupling Coefficients 

The time-dependent amplitude equations describe the evolution 
of the amplitudes of the excited modes. From the knowledge of 
the transient hydrodynamical tracks for a given model it is thus 
possible to evaluate the a priori unknown nonlinear coefficients of 
the amplitude equations. 

Of course, in order to extract reliable nonlinear coupling coeffi- 
cients it is necessary to integrate the models from a sufficient num- 
ber of initial conditions so that the 'trajectories' sample well the 
phase-phase (A ,Ai). Numerical simulations have show that mul- 
tiple DM solutions, a stable fixed point and a saddle point can co- 
exist in the model pulsations. From the analysis of the fixed-points 
of the amplitude equations one can conclude that it is necessary to 
include at least one quintic term into (see Section 2.). According to 
our experiments the mixed quintic coefficients (sj) are usually suf- 
ficient to give a good fit. The inclusion of the other quintic terms 
(rij) results only in a minor improvement of the fit. We have to 
note however that fits with the ru terms only are also capable of 
reproducing the main features of the amplitude evolution. 

In Table 1.3 we present the coefficients of the amplitude equations 
for a sequence of RR Lyrae models where a DM solution exists in 
the 6500.fr < T e ff < 6515/T temperature range. There are no 
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significant variations with T e ^ for most of the coefficients, but the 
growth-rate k has a decreasing trend. 



Table 1.3 Coupling coefficients in a sequence 
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6495 


0.009397 


0.03975 


-3.08 
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-44.20 


-2789. 


-729. 


6500 


0.009319 


0.03976 
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-20.29 
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-2788. 


-583. 


6505 


0.009235 


0.03967 


-3.07 


-20.35 


-10.13 


-43.94 


-2750. 


-500. 


6515 


0.009095 


0.03976 


-3.08 


-20.62 


-10.08 


-43.90 


-2833. 


-384. 


6525 


0.008896 


0.03971 


-3.07 


-20.75 


-9.99 


-43.72 


-2792. 


-360. 



We have fitted the coefficients with a quadratic function of tem- 
perature (e.g. sj = C + CiT e ff + C 2 T e 2 ^), and we have used the 
fit to map the solution of the amplitude equations as a function of 
temperature. Fig 1.7 shows the variation of the function 5F as the 
temperature is decreased. The displayed behavior is very similar 
to the one displayed in Fig 1.2, where only one of the parameters, 
the fundamental mode growth-rate was changed. If we fix all the 
parameters but decrease k according to the temperature fit, the 
basic tendencies on the figures remain the same i.e. for the given 
example the variation of the fundamental growth-rates plays the 
most important role. 

The Fourier representation of steady DM pulsations 

The primary information one can derive from observational data 
is the spectral content of the light variations, the amplitudes and 
periods presented in the Fourier spectra. (Examples of RR Lyrae 
frequency spectra can be found in Kovacs's review in this Volume.) 
In addition, because of the nonlinearity of pulsation, linear combi- 
nations of the primary frequencies also appears in the spectra. A 
good example for the existence of those terms is the analysis of the 
Cepheid TU Cassiopeia (Szabados 1993), where 29 frequencies have 
been identified. Even terms like 4fo+3fi or 5fo+2fi are inferred from 
the Fourier transform of the light-curve. In the RR Lyrae spectra 
typically only the fo+fi and fo fi terms can be found in addition to 
fo and fi. 
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Figure 1.7 Double-mode solutions of the amplitude equations fitted to a model se- 
quence. 

In order to relate the observed frequency spectra with those cal- 
culated from model pulsations, we have calculated synthetic data 
with the observed periodicities but with the same sampling and 
length as the hydrodynamical outputs (typically 1000 cycle). In 
Figure 1.8 we compare the spectra of TU Cas and of a MACHO 
DM star presented in Figure 2 of Kovacs in this Volume with typi- 
cal hydrodynamical results. We note that the amplitudes are given 
in bolometric magnitudes for the models and in 'V (TU Cas) and 
in instrumental red (MACHO RR Lyr). The agreement is satis- 
factory between the RR Lyrae models and observations. TU Cas 
has a large amplitude variation - it is the high nonlinearity that 
causes the large power at the harmonics of the fundamental mode. 
Our model has a much lower amplitude and consequently a weaker 
harmonic power, and it thus misses this behavior. 

When the harmonics of the modal frequencies can be observed, 
the corresponding Fourier parameters provide useful information. 
It has been observed that the Fourier parameters, y?2i and ip^i, 
of both modes in DM Cepheids are very similar to those of the 
corresponding modes in mono-periodic stars with similar periods 
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Figure 1.8 Fourier spectrum of the light variation: From top to bottom: TU Cas- 
siopeia, MACHO RRd variable, Cepheid model, RR Lyrae model. 



{e.g. Udalski et al. 1999, Beaulieu, private communication). Since 
the models relax very slowly on the separatrix connecting the F 
mode to the Oi through the DM solution, one can calculate the 
Fourier parameters along the separatrix. This test provides a clue 
how the Fourier phases depend on the modal amplitudes (Fig. 1.9). 
We have selected a model from the middle of the DM region for the 
plot, but we checked the Fourier parameters for the whole range 
and obtained very similar results. The Fourier parameters were 
calculated from the hydrodynamical model light-curves, and the 
amplitudes are in magnitudes. To indicate the position of the sep- 
aratrix, we also show the A\ vs. Aq plot of the models in the top 
right box. As expected the R21 amplitude ratios have a strong de- 
pendence on the corresponding amplitudes. For the fundamental 
mode this dependence is almost exactly linear. The plots of the ip 2 i 
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Fourier phases are very flat because there are no resonances. The 
difference between the </? 2 i value of the DM solution and that of the 
SM limit-cycle is less than 0.2. The y?3i parameters have the same 
tendencies but with a little bit stronger amplitude dependence for 
the F mode. 
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Figure 1.9 The dependence of the Fourier parameters on the fundamental mode 
amplitude along the separatrix. Filled hexagons: fundamental mode, open squares: 
first overtone. The top right box shows the models on the (Ao, Ai) phase-space. 



The Period ratios 

The period ratios offer an additional constraint for determin- 
ing the physical parameters of DM Cepheids and RR Lyrae stars 
(e.g. Kovacs in this Volume). These investigations make use of the 
linear periods and the latter are usually based on radiative pul- 
sation models. Thus it is very important to check how turbulent 
convection and nonlinear effects affect the period ratios. 
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Turbulent convection 

The dominant source of period shifts between the radiative and 
turbulent convective models is the convection induced change in 
the structure. In the radiative models there is a huge temperature 
gradient around the hydrogen partial ionization zone, because ra- 
diation is not sufficiently effective in transporting energy. In the 
turbulent convective models the convective flux carries some en- 
ergy and thus reduces the temperature gradient. This change in 
the stellar structure causes a decrease in the period ratios. The 
shift depends on the efficiency of the convection, and thus is larger 
for lower temperature models. Since the convective flux parameter 
(a c ) has not been satisfactorily calibrated so far, it is not possible to 
state definitive values for the shifts in the period ratios. Depending 
on the parameters of the models one finds (Pi/Po)tc ~ {P\/ Po)rad ~ 
-0.0005 to -0.002. Turbulent viscosity causes an additional shift of 
the periods, but that is small compared to the effect of the convec- 
tive flux. 

Nonlinearity 

The magnitude of the nonlinear period shifts has remained un- 
known due to the lack of nonlinear DM pulsation models. In order 
to isolate the period shifts from the larger period variations with 
T e ff we have used the time-frequency analysis of the model output 
together with the amplitude-equation formalism, instead of check- 
ing only the periods in the stationary DM oscillations. 

The phase that is calculated from the analytic signal (Eq. 1.37) 
is sufficiently smooth to yield the instantaneous frequencies by a 
numerical time derivative of y?(t). Hydrodynamical integrations 
with a set of different initial perturbations then can be used to fit 
the coefficients in the phase equation (Eq. 1.7). Since the nonlinear 
correction of the frequencies are small compared to Uj, the period 
ratios can be fitted directly with monomials 



(Pl/Pp)nl — (Pl/Po)lin 
{Pl/Po)lin 




(1.38) 



Fig. 1.10 and 1.11 shows examples of the hydrodynamical tracks 
that have been used to extract the nonlinear instantaneous period 
shift as a function of the amplitudes in the transient states. Fig 1.10 
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refers to an RR Lyrae model (M=O.77M , L=5OL ) and Fig 1.11 
to a Cepheid model (M=4.OM , L=11OOL ). Although this fit 
gives an interpolation formula for a specific model (with given mass, 
luminosity, and effective temperature) only, the fit is found to vary 
only weakly inside the DM region. 




Figure 1.10 Nonlinear period shift of an RR Lyrae model in the (Ao, Ai) phase-space. 

The 'equi-period-ratio' curves in the phase-space of the ampli- 
tudes are also shown in Fig 1.10 and 1.11. The relative changes of 
the nonlinear (TC) period ratios compared to the linear ones are 
indicated in the figures. For the Cepheid model the period ratio 
is always smaller than the linear one. In contrast, for RR Lyrae 
models the shift can be both positive and negative. These trends in 
the nonlinear shifts of the period-ratios have been found to persist 
in our model calculations throughout a wide range of model param- 
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eters. Since the vast majority of the RR Lyrae stars pulsate with 
higher amplitude in the first overtone than the fundamental mode, 
the observed period ratios are likely to be smaller than the linear 
values. 

Combining the two results we can conclude that the nonlinear 
turbulent convective period ratios are thus generally but not al- 
ways, smaller than the linear radiative ones, and that the relative 
difference is of the order of several tenths of a percent. 

Modal selection — Bifurcation Diagram 

The model calculations with different initial perturbations of the 
static model give the possible fixed point solutions of the star. The 
behavior of the individual models can be classified into five groups: 
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(1) first overtone only (O), 

(2) fundamental mode only (F), 

(3) either fundamental or first overtone (E), 

(4) DM only (D), and 

(5) either fundamental or DM (H). 

We have not encountered any model where the first overtone and 
the DM fixed points are simultaneously stable, but in principle such 
a situation is possible. For models with 'E' or 'H' characteristics 
there is hysteresis and one should examine a sequence along the 
possible evolutionary path of the star to predict its pulsation state. 
Sequences with different temperatures, but with a fixed luminosity, 
give good estimates for that purpose. 
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Figure 1.12 Modal selection on the (T e s, a„) plane for a Cepheid model. The symbols 
of the modal states are described at the top of the figure. 

In Fig. 1.12 the dependence of the modal selection on the eddy 
viscosity is shown for a Cepheid model sequence with SMC compo- 
sition. The Galactic Cepheids and the RR Lyrae models show very 
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similar behavior (for similar a parameters). The consistency of this 
similarity is somewhat astonishing in light of the serious discrepan- 
cies that one finds between the the observations and the calculations 
of low metallicity single-mode Cepheids (to be discussed elsewhere). 
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Figure 1.13 Modal selection of galactic Cepheid model sequences on the HR diagram. 
See Fig. 1.12 for the notation 



Turbulent viscosity has two dominant effects. The first one is to 
shift the transition region to higher temperatures (and to shorter 
periods) at higher a v values. The second one is to changes the 
characteristics of modal selection. For high viscosity parameters a 
DM-only solution exists sandwiched between the O and F state (we 
denote this transition by O-D-F). As the eddy viscosity decreased 
an 'either F or DM' (H) state appears at temperatures between 
the D and F solutions (O-D-H-F transition). For low a v values the 
characteristics of the transition changes to O-E-H-F i.e. the DM 
only region is changed to an 'either-or' one. 

For different stellar parameters (see Fig. 1.13) in addition to the 
transition types seen in Fig. 1.12 we could identify the following 
variations of modal behavior: O-E-F and O-D-H-E-F sequence. 

For the O-D-F type sequence DM behavior occurs independently 
of the direction of evolution. However, for the O-E-H-F case DM 
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pulsation can be observed only if the star evolves from the blue 
side of the instability strip to the red. In that case the star evolves 
to the DM state, when the Oi mode loses its stability, while on 
the opposite evolutionary path the star moves to the overtone di- 
rectly from the F mode at the blue edge of the fundamental mode 
instability strip through a transient DM state. 

The modal selection is illustrated in Figures 1.13 and 1.14 both 
for a Galactic Cepheid model and for an RR Lyrae model in an HR 
diagram . 

We have calculated the Cepheid models with a mass-luminosity 
relation (the masses of the models are indicated at the right end 
of the gray stripes). One notes a very interesting feature, namely 
the disappearance of the DM solution at higher luminosities, even 
where stable first overtone pulsation exists. This result is in good 
agreement with the observations by Udalski et al. (1999). 

Despite this very encouraging agreement we recall however that 
these bifurcation diagrams remain tentative until a full calibration 
of the a parameters has been made. 
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5. THE FUTURE 

Thanks to the incorporation of turbulent convection in the hy- 
drodynamical codes it is finally possible to model steady DM pulsa- 
tions both in Cepheids and RR Lyrae stars. Some problems persist 
though. The large-scale surveys have found second overtone to 
first overtone (O2/O1) DM pulsators at similar population as Oi/F 
beat Cepheids. The models should reproduce this constraint as 
well. With the recent parameter setting we could find some 2 /Oi 
models, but in a very narrow range of temperature and turbulence 
parameters. To model these DM variables satisfactorily, we have 
to explore the a parameter space in a wider range and perhaps 
improve the description of turbulent convection. 

This effort is important because a consistent modelling of DM 
oscillations in both Cepheids and RR Lyrae stars and with different 
chemical compositions is necessary for understanding the physical 
mechanisms inside the stars, on the one hand, and for enhancing our 
confidence that we can use these stars reliably as standard candles, 
on the other hand. 
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